<html>
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
    <title>leastsq</title>
  </head>
  <body bgcolor="#FFFFFF">
    <center>Scilab Function</center>
    <div align="right">Last update : 11/03/2005</div>
    <p>
      <b>leastsq</b> -  
    Solves non-linear least squares problems  
  </p>
    <h3>
      <font color="blue">Calling Sequence</font>
    </h3>
    <dl>
      <dd>
        <tt>[fopt,[xopt,[grdopt]]]=leastsq(fun, x0)</tt>
      </dd>
      <dd>
        <tt>[fopt,[xopt,[grdopt]]]=leastsq(fun, dfun, x0)</tt>
      </dd>
      <dd>
        <tt>[fopt,[xopt,[grdopt]]]=leastsq(fun, cstr, x0)</tt>
      </dd>
      <dd>
        <tt>[fopt,[xopt,[grdopt]]]=leastsq(fun, dfun, cstr, x0)</tt>
      </dd>
      <dd>
        <tt>[fopt,[xopt,[grdopt]]]=leastsq(fun, dfun, cstr, x0, algo)</tt>
      </dd>
      <dd>
        <tt>[fopt,[xopt,[grdopt]]]=leastsq([imp], fun [,dfun] [,cstr],x0 [,algo],[df0,[mem]],[stop])</tt>
      </dd>
    </dl>
    <h3>
      <font color="blue">Parameters</font>
    </h3>
    <ul>
      <li>
        <tt>
          <b>fopt</b>
        </tt>: value of the function <em>f(x)=||fun(x)||^2</em> at <tt>
          <b>xopt</b>
        </tt>
      </li>
      <li>
        <tt>
          <b>xopt</b>
        </tt>: best value of <tt>
          <b>x</b>
        </tt> found to minimize <em>||fun(x)||^2</em>
      </li>
      <li>
        <tt>
          <b>grdopt</b>
        </tt>: gradient of <em>f</em> at <tt>
          <b>xopt</b>
        </tt>
      </li>
      <li>
        <tt>
          <b>fun</b>
        </tt>:  a scilab function or a list defining a function from <em>R^n</em> to <em>R^m</em> 
                 (see more details in DESCRIPTION).</li>
      <li>
        <tt>
          <b>x0</b>
        </tt>: real vector (initial guess of the variable to be minimized).</li>
      <li>
        <tt>
          <b>dfun</b>
        </tt>:  a scilab function or a string defining the Jacobian matrix of <tt>
          <b>fun</b>
        </tt>
                 (see more details in DESCRIPTION).</li>
      <li>
        <tt>
          <b>cstr</b>
        </tt>: bound constraints on <tt>
          <b>x</b>
        </tt>. They must be introduced by the string keyword <tt>'b'</tt>
                followed by the lower bound <tt>binf</tt> then by the upper bound <tt>bsup</tt> (so 
                <tt>
          <b>cstr</b>
        </tt> appears as <tt>
          <b>'b',binf,bsup</b>
        </tt> in the calling sequence). Those bounds
                are real vectors with same dimension than <tt>
          <b>x0</b>
        </tt> (-%inf and +%inf may be used for
                dimension which are unrestricted).</li>
      <li>
        <tt>
          <b>algo</b>
        </tt>: a string with possible values: <tt>'qn'</tt> or <tt>'gc'</tt> or <tt>'nd'</tt>. 
                These strings stand for quasi-Newton (default),  conjugate gradient or non-differentiable 
                respectively. Note that <tt>'nd'</tt> does not accept bounds on <tt>
          <b>x</b>
        </tt>.</li>
      <li>
        <tt>
          <b>imp</b>
        </tt>: scalar argument used to set the trace mode. <tt>
          <b>imp=0</b>
        </tt> nothing (except errors) 
                is reported, <tt>
          <b>imp=1</b>
        </tt> initial and final reports, <tt>
          <b>imp=2</b>
        </tt> adds a 
                report per iteration, <tt>
          <b>imp&gt;2</b>
        </tt> add reports on linear search. Warning, most
                of these reports are written on the Scilab standard output.</li>
      <li>
        <tt>
          <b>df0</b>
        </tt>: real scalar. Guessed decreasing of <em>||fun||^2</em> at first iteration. 
                (<tt>df0=1</tt> is the default value).</li>
      <li>
        <tt>
          <b>mem</b>
        </tt>: integer, number of variables used to approximate the Hessean (second derivatives)
                of <em>f</em> when <tt>
          <b>algo</b>
        </tt>
        <tt>='qn'</tt>. Default value is around 6.</li>
      <li>
        <tt>
          <b>stop</b>
        </tt>:  sequence of optional parameters controlling the  convergence of the algorithm. They 
                 are introduced by the keyword <tt>'ar'</tt>, the sequence being of the form
                 <tt>
          <b>'ar',nap, [iter [,epsg   [,epsf [,epsx]]]]</b>
        </tt>
        <ul>
          <li>
            <tt>
              <b>nap</b>
            </tt>: maximum number of calls to <tt>
              <b>fun</b>
            </tt> allowed.</li>
          <li>
            <tt>
              <b>iter</b>
            </tt>: maximum number of iterations allowed.</li>
          <li>
            <tt>
              <b>epsg</b>
            </tt>: threshold on gradient norm.</li>
          <li>
            <tt>
              <b>epsf</b>
            </tt>: threshold controlling decreasing of <tt>
              <b>f</b>
            </tt>
          </li>
          <li>
            <tt>
              <b>epsx</b>
            </tt>: threshold controlling variation of <tt>
              <b>x</b>
            </tt>. This vector (possibly matrix) 
                      of same size as <tt>
              <b>x0</b>
            </tt> can be used to scale <tt>
              <b>x</b>
            </tt>.</li>
        </ul>
      </li>
    </ul>
    <h3>
      <font color="blue">Description</font>
    </h3>
    <p>
      <em>fun</em> being a function from <em>R^n</em> to <em>R^m</em> this routine tries to minimize w.r.t.
     <em>x</em>, the function:
   </p>
    <pre>
                           _m_     
                      2    \      2
     f(x) = ||fun(x)||  =  /   fun (x)
                           ---    i
                           i=1     
          </pre> which is the sum of the squares of the components of <em>fun</em>. Bound constraints may be
     imposed on <tt>
      <b>x</b>
    </tt>.<h3>
      <font color="blue">How to provide fun and dfun</font>
    </h3>
    <dl>
      <p>
        <tt>
          <b>fun</b>
        </tt> can be either a usual scilab function (case 1) or a fortran or a C routine linked
            to scilab (case 2). For most problems the definition of <em>fun</em> will need 
            supplementary parameters and this can be done in both cases.</p>
      <dd>
        <li>
          <b>
            <font color="maroon">case 1:</font>
          </b> when <tt>
            <b>fun</b>
          </tt> is a Scilab function, its calling sequence must be:
          <tt>y=fun(x [,opt_par1,opt_par2,...])</tt>. When <tt>
            <b>fun</b>
          </tt> needs optional parameters it
           must appear as <tt>list(fun,opt_par1,opt_par2,...)</tt> in the calling sequence of
           <tt>
            <b>leastsq</b>
          </tt>.</li>
        <li>
          <b>
            <font color="maroon">case 2:</font>
          </b> when <tt>
            <b>fun</b>
          </tt> is defined by a Fortran or C routine it must appear as
                <tt>list(fun_name,m [,opt_par1,opt_par2,...])</tt> in the calling
                sequence of <tt>
            <b>leastsq</b>
          </tt>, <tt>fun_name</tt> (a string) being the name 
                of the routine which must be linked to Scilab (see <a href="../utilities/link.htm">
            <tt>
              <b>link</b>
            </tt>
          </a>). The generic 
                calling sequences for this routine are:<pre>
In Fortran:    subroutine fun(m, n, x, params, y)
               integer m,n
               double precision x(n), params(*), y(m)

In C:          void fun(int *m, int *n, double *x, double *params, double *y)
         </pre>where <tt>
            <b>n</b>
          </tt> is the dimension of vector <tt>
            <b>x</b>
          </tt>, <tt>
            <b>m</b>
          </tt> the dimension of 
            vector <tt>
            <b>y</b>
          </tt> (which must store the evaluation of <em>fun</em> 
            at <em>x</em>) and <tt>
            <b>params</b>
          </tt> is a vector which contains the optional parameters 
            <em>opt_par1, opt_par2, ...</em> (each parameter may be a vector, for instance if <em>opt_par1</em> 
            has 3 components, the description of <em>opt_par2</em> begin from <tt>params(4)</tt> (fortran case),
            and from  <tt>params[3]</tt> (C case), etc... Note that even if <tt>
            <b>fun</b>
          </tt> doesn't need
            supplementary parameters you must anyway write the fortran code with a <tt>params</tt> argument
            (which is then unused in the subroutine core).</li>
      </dd>
      <p>In many cases it is adviced to provide the Jacobian matrix <tt>
          <b>dfun</b>
        </tt> (<em>dfun(i,j)=dfi/dxj</em>) 
           to the optimizer (which uses a finite difference approximation otherwise) and as for <tt>
          <b>fun</b>
        </tt> it may 
           be given as a usual scilab function or as a  fortran or a C routine linked to scilab.</p>
      <dd>
        <li>
          <b>
            <font color="maroon">case 1:</font>
          </b> when <tt>
            <b>dfun</b>
          </tt> is a scilab function, its calling sequence must be:
          <tt>y=dfun(x [, optional parameters])</tt> (notes that even if <tt>
            <b>dfun</b>
          </tt> needs optional 
              parameters it must appear simply as <tt>dfun</tt> in the calling sequence of
           <tt>
            <b>leastsq</b>
          </tt>).</li>
        <li>
          <b>
            <font color="maroon">case 2:</font>
          </b> when <tt>
            <b>dfun</b>
          </tt> is defined by a Fortran or C routine it must appear as
                <tt>dfun_name</tt> (a string) in the calling sequence of <tt>
            <b>leastsq</b>
          </tt> 
                (<tt>dfun_name</tt> being the name of the routine which must be linked to Scilab). 
                The calling sequences for this routine are nearly the same than for <tt>
            <b>fun</b>
          </tt>:<pre>
In Fortran:    subroutine dfun(m, n, x, params, y)
               integer m,n
               double precision x(n), params(*), y(m,n)

In C:          void fun(int *m, int *n, double *x, double *params, double *y)
         </pre>in the C case <em>dfun(i,j)=dfi/dxj</em> must be stored in <tt>y[m*(j-1)+i-1]</tt>.</li>
      </dd>
    </dl>
    <h3>
      <font color="blue">Remarks</font>
    </h3>
    <dl>
      <p> Like <a href="datafit.htm">
          <tt>
            <b>datafit</b>
          </tt>
        </a>, <tt>
          <b>leastsq</b>
        </tt> is a front end onto the <a href="optim.htm">
          <tt>
            <b>optim</b>
          </tt>
        </a>
       function. If you want to try the Levenberg-Marquard method instead, use <a href="lsqrsolve.htm">
          <tt>
            <b>lsqrsolve</b>
          </tt>
        </a>.
    </p>
      <p> A least squares problem may be solved directly with the <a href="optim.htm">
          <tt>
            <b>optim</b>
          </tt>
        </a> function ;
        in this case the function <a href="NDcost.htm">
          <tt>
            <b>NDcost</b>
          </tt>
        </a> may be useful to compute the derivatives
        (see the <a href="NDcost.htm">
          <tt>
            <b>NDcost</b>
          </tt>
        </a> help page which provides a simple example for parameters 
        identification of a differential equation).
    </p>
    </dl>
    <h3>
      <font color="blue">Examples</font>
    </h3>
    <pre>
// We will show different calling possibilities of leastsq on one (trivial) example
// which is non linear but doesn't really need to be solved with leastsq (applying
// log linearizes the model and the problem may be solved with linear algebra). 
// In this example we look for the 2 parameters x(1) and x(2) of a simple
// exponential decay model (x(1) being the unknow initial value and x(2) the
// decay constant): 
function y = yth(t, x)
   y  = x(1)*exp(-x(2)*t) 
endfunction  

// we have the m measures (ti, yi):
m = 10;
tm = [0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5]';
ym = [0.79, 0.59, 0.47, 0.36, 0.29, 0.23, 0.17, 0.15, 0.12, 0.08]';
wm = ones(m,1); // measure weights (here all equal to 1...)

// and we want to find the parameters x such that the model fits the given 
// datas in the least square sense:
// 
//  minimize  f(x) = sum_i  wm(i)^2 ( yth(tm(i),x) - ym(i) )^2   

// initial parameters guess
x0 = [1.5 ; 0.8];

// in the first examples, we define the function fun and dfun 
// in scilab language
function e = myfun(x, tm, ym, wm)
   e = wm.*( yth(tm, x) - ym )
endfunction

function g = mydfun(x, tm, ym, wm)
   v = wm.*exp(-x(2)*tm)
   g = [v , -x(1)*tm.*v]
endfunction

// now we could call leastsq:

// 1- the simplest call
[f,xopt, gropt] = leastsq(list(myfun,tm,ym,wm),x0)

// 2- we provide the Jacobian
[f,xopt, gropt] = leastsq(list(myfun,tm,ym,wm),mydfun,x0)

// a small graphic (before showing other calling features)
tt = linspace(0,1.1*max(tm),100)';
yy = yth(tt, xopt);
xbasc()
plot2d(tm, ym, style=-2)
plot2d(tt, yy, style = 2)
legend(["measure points", "fitted curve"]);
xtitle("a simple fit with leastsq")

// 3- how to get some informations (we use imp=1)
[f,xopt, gropt] = leastsq(1,list(myfun,tm,ym,wm),mydfun,x0)

// 4- using the conjugate gradient (instead of quasi Newton)
[f,xopt, gropt] = leastsq(1,list(myfun,tm,ym,wm),mydfun,x0,"gc")

// 5- how to provide bound constraints (not useful here !)
xinf = [-%inf,-%inf]; xsup = [%inf, %inf];
[f,xopt, gropt] = leastsq(list(myfun,tm,ym,wm),"b",xinf,xsup,x0) // without Jacobian
[f,xopt, gropt] = leastsq(list(myfun,tm,ym,wm),mydfun,"b",xinf,xsup,x0) // with Jacobian 

// 6- playing with some stopping parameters of the algorithm
//    (allows only 40 function calls, 8 iterations and set epsg=0.01, epsf=0.1)
[f,xopt, gropt] = leastsq(1,list(myfun,tm,ym,wm),mydfun,x0,"ar",40,8,0.01,0.1)


// 7 and 8: now we want to define fun and dfun in fortran then in C
//          Note that the "compile and link to scilab" method used here
//          is believed to be OS independant (but there are some requirements, 
//          in particular you need a C and a fortran compiler, and they must 
//          be compatible with the ones used to build your scilab binary).
 
// 7- fun and dfun in fortran

// 7-1/ Let 's Scilab write the fortran code (in the TMPDIR directory):
f_code = [ ...
"      subroutine myfun(m,n,x,param,f)"
"*     param(i) = tm(i), param(m+i) = ym(i), param(2m+i) = wm(i)"
"      implicit none"
"      integer n,m"
"      double precision x(n), param(*), f(m)"
"      integer i"
"      do i = 1,m"
"         f(i) = param(2*m+i)*( x(1)*exp(-x(2)*param(i)) - param(m+i) )"
"      enddo"
"      end ! subroutine fun"
""
"      subroutine mydfun(m,n,x,param,df)"
"*     param(i) = tm(i), param(m+i) = ym(i), param(2m+i) = wm(i)"
"      implicit none"
"      integer n,m"
"      double precision x(n), param(*), df(m,n)"
"      integer i"
"      do i = 1,m"
"         df(i,1) =  param(2*m+i)*exp(-x(2)*param(i))"
"         df(i,2) = -x(1)*param(i)*df(i,1)"
"      enddo"
"      end ! subroutine dfun"];
mputl(f_code,TMPDIR+'/myfun.f')

// 7-2/ compiles it. You need a fortran compiler !
names = ["myfun" "mydfun"]
flibname = ilib_for_link(names,"myfun.o",[],"f",TMPDIR+"/Makefile");

// 7-3/ link it to scilab (see link help page)
link(flibname,names,"f") 

// 7-4/ ready for the leastsq call: be carreful don't forget to
//      give the dimension m after the routine name !
[f,xopt, gropt] = leastsq(list("myfun",m,tm,ym,wm),x0)  // without Jacobian
[f,xopt, gropt] = leastsq(list("myfun",m,tm,ym,wm),"mydfun",x0) // with Jacobian


// 8- last example: fun and dfun in C

// 8-1/ Let 's Scilab write the C code (in the TMPDIR directory):
c_code = [...
"#include &lt;math.h&gt;"
"void myfunc(int *m,int *n, double *x, double *param, double *f)"
"{"
"  /*  param[i] = tm[i], param[m+i] = ym[i], param[2m+i] = wm[i] */"
"  int i;"
"  for ( i = 0 ; i &lt; *m ; i++ )"
"    f[i] = param[2*(*m)+i]*( x[0]*exp(-x[1]*param[i]) - param[(*m)+i] );"
"  return;"
"}"
""
"void mydfunc(int *m,int *n, double *x, double *param, double *df)"
"{"
"  /*  param[i] = tm[i], param[m+i] = ym[i], param[2m+i] = wm[i] */"
"  int i;"
"  for ( i = 0 ; i &lt; *m ; i++ )"
"    {"
"      df[i] = param[2*(*m)+i]*exp(-x[1]*param[i]);"
"      df[i+(*m)] = -x[0]*param[i]*df[i];"
"    }"
"  return;"
"}"];
mputl(c_code,TMPDIR+'/myfunc.c')

// 8-2/ compiles it. You need a C compiler !
names = ["myfunc" "mydfunc"]
clibname = ilib_for_link(names,"myfunc.o",[],"c",TMPDIR+"/Makefile");

// 8-3/ link it to scilab (see link help page)
link(clibname,names,"c") 

// 8-4/ ready for the leastsq call
[f,xopt, gropt] = leastsq(list("myfunc",m,tm,ym,wm),"mydfunc",x0)
   
  </pre>
    <h3>
      <font color="blue">See Also</font>
    </h3>
    <p>
      <a href="lsqrsolve.htm">
        <tt>
          <b>lsqrsolve</b>
        </tt>
      </a>,&nbsp;&nbsp;<a href="optim.htm">
        <tt>
          <b>optim</b>
        </tt>
      </a>,&nbsp;&nbsp;<a href="NDcost.htm">
        <tt>
          <b>NDcost</b>
        </tt>
      </a>,&nbsp;&nbsp;<a href="datafit.htm">
        <tt>
          <b>datafit</b>
        </tt>
      </a>,&nbsp;&nbsp;<a href="../programming/external.htm">
        <tt>
          <b>external</b>
        </tt>
      </a>,&nbsp;&nbsp;<a href="quapro.htm">
        <tt>
          <b>quapro</b>
        </tt>
      </a>,&nbsp;&nbsp;<a href="linpro.htm">
        <tt>
          <b>linpro</b>
        </tt>
      </a>,&nbsp;&nbsp;</p>
  </body>
</html>
